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Abstract. We derive recursions for the probability distribution of random 
sums by computer algebra. Unlike the well-known Panjer-type recursions, 
they are of finite order and thus allow for computation in linear time. This 
efficiency is bought by the assumption that the probability generating function 
of the claim size be algebraic. The probability generating function of the claim 
number is supposed to be from the rather general class of D-finite functions. 
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1. Introduction 

Random sums 

L = Xi + • • • + Xj\i , 

play a prominent role in risk theory. We refer to the X^, which are independent 
copies of a discrete random variable X, as claims, and to N , which is indepen- 
dent of the Xi, as claim number. A lot of research has been devoted to recursive 
calculation of the distribution of L. The classical Panjer recursion and its numer- 
ous extensions |18l [71 \T'd\ provide infinite order linear recursions for this problem 
for various claim number distributions. Recently, Hipp [5] has found that finite 
order recursions can be obtained for phase-type claim size distributions, with obvi- 
ous advantages concerning computation time. The generalized discrete phase-type 
distributions are the distributions of hitting times in a finite-state discrete-time 
Markov chain. They serve as a very flexible class of severity distributions with sev- 
eral useful properties [T^ [5] • Simple examples are the geometric and the negative 
binomial distribution (with integral parameter a). Another approach by DePril [3] 
also leads to finite order recursions for some claim size distributions, for example for 
piecewise constant or piecewise linear claim size distributions. But the recursion by 
DePril is not of finite order for the broad class of algebraic probability generating 
functions which we cover in this article. 

We present an alternative method to obtain recursions of finite order satisfied 
by the distribution P[L = n] of L. Our assumption on the claim number N is 
that its probability generating function (pgf) be a Z?-finite function. This broad 
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class of functions is characterized by linear differential equations with polynomial 
coefficients. The pgf of the claims is assumed to be algebraic. Then the theory 
of _D-finite functions ensures the existence of a finite order linear recursion with 
polynomial coefficients for the distribution of L. The recursion can be determined 
by computer algebra, which avoids tedious hand calculations. 

Section [2] collects the parts of the theory of _D-finite functions that we want 
to use, together with their algorithmic realization in computer algebra systems. 
The latter is applied in Section [3] to some concrete examples. Although we focus 
on computational efficiency, we do care about numerical stability of the proposed 
recursions. Therefore, in Section |4] we apply the stability theory of finite order 
recursions to our examples. 

2. _D-FINITE FUNCTIONS 

We will assume throughout that the probability generating (pgf) function 
of the claim number N is of the following kind. 

Definition 2.1. Let f{z) be a function that is analytic at zero. Then f{z) is called 
D- finite if it satisfies a linear differential equation 

Qo{z)f{z) + Qi{z)f\z) + ■■■ + Qd{z)f'^^\z) = (2.1) 

with polynomial coefficients Qo{z), . . . , Qd{z), not all identically zero. 

Most discrete distributions that are used in practice have _D-finite pgfs. See, 
e.g., the comprehensive list of hypergeometric distributions in Johnson, Kemp, 
Kotz [TD] . Many properties of D-finite functions follow from the classical theory of 
ordinary differential equations [9] ; see Stanley [16l [17] for an introduction from a 
combinatorial viewpoint including a proof of the following result. 

Theorem 2.2. 

(i) An analytic function f{z) — X)n>o "^n^" D-finite if and only if its coeffi- 
cient sequence {an)n>o satisfies a finite order linear recursion with polyno- 
mial coefficients. 

(ii) Ttie sum and tfie product of two D-finite functions are D-finite. The com- 
position f{g{z)) of a D-finite function f{z) and an algebraic function g{z) 
is D-finite. 

Part (i) says that a differential equation of the form (|2.ip always translates into 
a recursion 

Ro{n)an + i?i(n)a„+i H h Re{n)an+e =0, n > 0, 

with polynomial coefficients Rk(n) for the power series coefficients (a„)„>o of /(z), 
and vice versa. Therefore, the distributions from the Panjer class are simple ex- 
amples of distributions with _D-finite pgf. More examples of such distributions can 
be found, e.g., in Panjer and Willmot [13]. Hesselager [7] has found a recursion for 
the distribution of L for claim number distributions with arbitrary £)-finite pgf. It 
is valid for any claim size distribution, but is of infinite order in general. We, on 
the other hand, aim at finite order recursions for the distribution of L. Part (ii) 
of Theorem 12.21 is central for our approach. Recall that an algebraic function g{z) 
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satisfies P{z,g{z)) = for some non-trivial bivariate polynomial P. We explicitly 
note the result that part (ii) of Theorem 12.21 implies in our situation. 

Corollary 2.3. // the pgf lpn{z) of the claim numbers is D-finite and the pgf ipx{z) 
of the severity distribution is algebraic, then 



n>0 

is D-finite, and its coefficients a„ satisfy a linear recurrence of finite order with 
polynomial coefficients. 

Examples of admissible severity distributions include those whose pgfs are poly- 
nomials or rational functions, in particular, the phase- type distributions mentioned 
in the introduction. Furthermore, consider the negative binomial distribution 
NBin(a,p) with a > and p € (0, 1), i.e., 



is algebraic if a is a rational number. Other examples of distributions with algebraic 
pgf are the binomial distribution, the discrete Mittag-Leffler distribution, and its 
generalization, the discrete Linnik distribution [10'. Again, parameters appearing 
in the exponent have to be constrained to the rational numbers. The distribution 
of the number of games lost by the ruined gambler in the classical gambler's ruin 
problem fW has an algebraic pgf, too. 

All operations described in Theorem 12.21 are constructive and have been imple- 
mented in computer algebra packages. We have done the computations below with 
Mallinger's Mathematica package GeneratingFunctions [IT. Another possible choice 
would have been Salvy and Zimmcrmann's Maple package gfun [15 . The software 
requires an algebraic equation for the claim size pgf and a differential equation for 
the claim number pgf. We assume that our claim size pgf is an explicit algebraic 
(or even rational) function, so the first of these two equations is obvious in our 
examples. As for the second one, suppose we have a closed form expression for the 
D-finite claim number pgf. A differential equation for it can be built up step by 
step, by starting from obvious differential equations (for the exponential function, 
say) and using commands that realize the closure properties from Theorem 12.21 
Finally, the package allows to convert the differential equation that we have thus 
found for fL{z) into the desired finite order recursion for its coefficients. 

The package GeneratingFunctions can deal with undetermined parameters. For 
instance, the differential equation (az + b)f'{z) — cz'^ f{z) = with parameters 
a, 6, c would be a valid input, and we could, e.g., compute a recurrence relation for 
the power series coefRcients of f{z). When inputting algebraic relations, however, 
exponents must be fixed: We can specify /(z)^ = {az + 6)'^, e.g., but not f{z)°' — 
{az + bY with undetermined a. 





Its pgf 




3. Examples 



We illustrate our approach by three examples. The intermediate differential 
equations obtained during the process are not displayed. First we consider N ^ 
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NBin(a,p) and X ~ NBin(/?, g) with p,q £ (0, 1). Then the pgf of the aggregate 
loss is given by 

Our goal is to find a differential equation for fLiz) and thence a recurrence relation 
for a„ = P[L = n]. As a byproduct of the stability analysis in Section 21 we will 
obtain the asymptotics of a„ for general a and /3. To compute a recurrence for 
a„ by computer algebra, however, these parameters have to be concrete rational 
numbers (see the last paragraph of Section [J) . We choose en ^ ^ and /? = ^ . Then 
the function Lpf^[z) satisfies the algebraic equation 

(1 - (1 ~p)z)ipN{zf 

from which the command AlgebraicEquationToDifferentialEquation computes a dif- 
ferential equation for lpn{z). From the latter and the algebraic equation 

{l-{\-q)z)^x{zf 

of i^x{z), the command Algebra icCom pose (cf. Corollarv l2.3p derives the differential 
equation 

~b{p-lf{q-lfqf{z) + 

10((7 - 1)2(1 - z - qz){l& -q + Zpq- 3p^q + p^q - 16z + 16qz)f'{z) + 

36{q - 1)(1 -z + qz)^{8 -3q + 9pq - 9p^q + ip^q - 8z + %qz)f"{z)+ ^'^''^^ 

72(1 - z + qzf{l - q3pq ~ 3p^q +p^q~ z + qz)f"'{z) = 

for ^pl{z) = ipNi^xiz))- Finally, this equation is transformed into the recurrence 

8n{l + 3n){2 + 3n){q - l)^a„+ 

{q - l)-'^(320 + 896n + 86471^ + 28871^ - 5q - iGnq - lOSn'^q - 72n^q + 15pq + 138npq + 32An^pq+ 
216n^pq - Ihp^q - I38np^q - 32An^p^q - 216nV9 + ^P^q + 46np3g + 108nV9 + 72nV9)an+i + 
(4 + 2n){q - 1)^(512 + 648n + 216^^ - ll3q - 21671^ - 113g - 2l&nq - IQSn^q^ 
339pg + 648?7re + 32'in^pq - 339p^q - M?>np^q - 32An^p^q + 113p^g + 2l&np^q + l{)%n^p^q)an+2+ 
36(2 + n)(3 + n){q - 1)(16 + 8n - 9g - Qnq + 27pq + I'&npq - 27p^q ~ ISnp'^q + 9p^q + 6np^q)a„+3 
= 72(2 + n){3 + n){A + n){-l + q- 3pq + 3p^q - p^g)a„+4 



(3.3) 

for the a„ by DifferentialEquationToRecurrenceEquation. Using this recurrence, the 
probability a„ can be computed with 0{n) operations. We will show in Section [4] 
that the computation is numerically stable. 

In our second example we suppose that N ^ Poisson(A) and A ~ GlG{tp, x, 
the generalized inverse Gaussian distribution with parameters 9 and i/'jX > 0. In 
this case we have 

^n(z) - -{^ + 2- 2z)-'/^ ■ Ke{^x{^ + 2-2z)), (3.4) 

where Kg{z) is a modified Bessel function of the second kind [1, p. 374]. Provided 
that is a rational number, the second factor is an algebraic function, hence D- 
finite. (Below we will iix6 — |.) The Bessel function Kg{z) is D-finite for any 9, by 
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virtue of its classical second order differential equation. Therefore, by Theorem l2.21 
our (pn{z) is indeed a Z?- finite function. 

As for the severities, we take them to be shifted geometrically distributed: X ^ 
Geo(l, q) with q e (0, 1), so that 

'Px(z) = ^. 

1 - (1 - q)z 

Once again we want to find a differential equation, and thence a recurrence for the 
power series coefficients, for the function 

(Pl{z) = (Pn{(Px{z)) =: const • f{ipx{z)) ■ Kg{g{(px{z))), (3.5) 

where f{z) and g{z) are algebraic functions defined according to (|3.4p . The com- 
mand AlgebraicEquationToDifferentialEquation computes a differential equation for 
f{ipx{z)) from the algebraic equation 

f{^x{z)r'^'^i^ + 2- . 

1 - (1 - q)z 

As mentioned above, this works for any rational 6; to perform the calculation 
step, we have to fix its value though, say 9 — ^. A differential equation for 
K2/3{g{ipx{z))) can be found with AlgebraicCompose. It takes as input the dif- 
ferential equation of 1^2/3(2;) and the algebraic equation 

Now that we have a differential equation for each of the two (non-constant) factors 
in (|3.5p . the command DECauchy computes a differential equation for (pN{(px{z)), 
which is transformed into a recursion for its power series coefficients a„ by DifFer- 
entialEquationToRecurrenceEquation. The recurrence we find is 

^3/ 



3n(l -I- n)(-l + qy{~2 - ^j; + i;q) ■ 



2(1 + n)(-l + qf{-l8 - I2n - 9-0 - 6nV + q + 3nq + Qijjq + Gmljq) ■ a„ 



' {U4:+U4n + 36n'^ + 72ip + 72nip + 18n'^ip-188q-202nq-54n^q-lUipq-lUnipq 
- 36n^ipq + Uq^ - 3xq^ + 58nq^ + 18n^q^ + 72^/jq^ + 72n^/jq^ + ISn^g^) ' an+2 
-t- 2(3 + n)(-30 - 12n - Ibip - 6ml; + 19q + 9nq + 15iljq + Gmpq) ■ a„+3 

+ 3(3 + n)(4 + n)(2-|-V') •a„+4 =0. (3.6) 



Summing up, if N has a mixed Poisson distribution Poisson(A) with A ~ GlG{ip, x, f 
and X ^ Geo(l, q), then we can compute the total loss probabilities with 0(n) op- 
erations by the recurrence p.6p . 

The third example we consider is ~ Poisson(A) and X ~ NBin(i,p) with 
p S (0, 1). To determine a differential equation for the pgf 



we use again the command AlgebraicCompose. Its input are the algebraic equation 

{l-il^p)z)ipxiz)^ =P 

and the differential equation 

'p'n{z) = ^^fNiz). 



P 



1/2 
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Algebra icCom pose then finds the differential equation 

A2p(l -p)Vl(z) + 6(1 - (1 -p)z)Vl W - 4(1 - (1 -p)z)yi{z) = 0. (3.7) 

Using DifFerentialEquationToRecurrenceEquation, we obtain the recurrence 

2n{2n + - pf an - (1 -pf{~X^p+ \2n^ + 24n + 12)a„+i 

+ 6(n + 2)(2n + 3)(1 - p)a„+2 = 4(n + 2){n + 3)a„+3 (3.8) 
for the probabihties a„. 

4. Numerical Stability and Asymptotics 

The computation of a sequence by a linear recurrence relation is numerically 
stable if the sequence grows at least as fast as any other solution of the recur- 
rence [201 113] . This is intuitively clear, since rounding errors will always add a 
portion of each member of a fundamental system of the recurrence to the solution 
we are computing. Asymptotically dominant solutions will therefore wipe out sub- 
ordinate solutions in the long run. In this section we show how to apply methods 
from asymptotic analysis to the examples from Section [3l The growth of the coef- 
ficients a„ depends on the location and nature of the singularity of the generating 
function (flLiz) that is closest to the origin To assess the growth of the other 
solutions of our recurrences, we have to analyze the dominating singularity of the 
differential equation for tpL(z). If it is regular, then a fundamental system can in 
principle be determined by Frobenius' method. Flajolet and Odlyzko's singularity 
analysis [3 [6j allows to obtain the growth rate of the power series coefficients of 
these solutions. Then, hopefully, we can read off that the solution we are interested 
in dominates the other ones. This fairly general method works in the first two 
examples from Section [3l 

Proposition 1. Let N ^ 'NBm{a,p) and X ^ NBin(/?, q) with p,q E (0, 1). Then 
the probabilities an = ^[L = n] satisfy 

an ~ Czj-"n"-i (4.1) 

as n oo, where 

l-g(l-p)l/^ {pq)^(l-p)^/f^ 
Z\ — and C — 



1-q r(a)/3"(l - g(l -p)i//3)"- 

Proof. The dominating singularity of (Pl{z) is located at z = zi. Moving the 
singularity to z = 1 and putting c := 1 — q{l — pY^^ , we find 

/ 

Vl{zzi) = 



From the expansion 



" ' q V f c 



1 — cz/ \1 — cy \ 1 — c 
1 / Pc 



(l_^) + 0((l-z)^) 
I — p \ 1 — c 
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we thus obtain 



- -(l_z)-" + 0((l-z)-"+i) (4.2) 



/3"(l-g(l-p)i//3)" 

as z tends to 1. The n-th power series coefficient of (1 — z)~" asymptotically equals 
n"^^/r(a) [6, Chapter VI]. The coefficients of the error term in (|4.2p are of smaller 
order [6l Chapter VI], whence the desired result. We note that an asymptotic 
expansion to arbitrary order can be obtained in the same way. □ 

To be able to derive a concrete recursion we had to assign values to the param- 
eters OL and (5 in Section [31 

Proposition 2. Lei TV NBin(i,p) and X NBin(i, g) with p,q e (0,1). Then 
the computation of the a„ by the recursion (|3.3|) is numerically stable. 

Proof. We show, by Frobenius' theory of power series solutions of differential equa- 
tions [9], that no solution of (|3.3p grows faster than (|4.ip (with a — ^ and f] — 
Instead of working directly with (|3.2p . we transform (|3.3p into a differential equa- 
tion. This is necessary because new solutions might creep in when passing from a 
differential equation to a recurrence or vice versa. The new differential equation, 
called £ in what follows, is solved by all generating functions of solutions of (|3.3p . 
It is of order four and has the same leading coefficient as ()3.2p . By equating this 
coefficient to zero, we find the dominating singularity zi ~ {1 — q{l — p)^)/{l — q). 
The indicial polynomial of £ is found by plugging in a generalized power series 
(z — Zl)"" X]fe°=o ^k{z — Zl)*"' with undetermined w and Ak and equating the coef- 
ficient of the lowest power of zi — z to zero. The lowest power turns out to be 
(z — zi)'^~^ with the coefficient {w — 2){w — l)w{l + 2w) times a constant. The 
degree of this indicial polynomial equals the order of the differential equation, hence 
Zl is a regular singularity of £. The roots of the indicial polynomial are the possible 
values of the exponent w in the generalized power series solution. The root w = —-^ 
leads, by singularity analysis, to a solution of £ whose coefficients grow like (|4.ip 
(with a — ^ and (3 — i). The other solutions of the fundamental system that 
Frobenius' method yields have either no singularity at zi or a logarithmic singu- 
larity at Zl. The latter type occurs because some roots of the indicial polynomial 
differ by integers, and the coefficients of the corresponding solutions of £ grow like 
1 /n times a power of log n fG'i Chapter VI] . □ 

The second example from Section |3] can be treated analogously: 

Proposition 3. Let N ~ Poisson(A) and A ^ GIG('0, Xj^)i where V', X and 9 
are positive. Furthermore assume that X ^ Geo(l,<7) with q G (0,1). Then the 
probabilities a„ = P[L = n] satisfy 

a„ ~ Cx-'/^D-%2nf-hi'' (4.3) 

as n oo, where 

1 , ^ (2 + V)(2 + v^(l-g)) 

^ = Zl = ■ — 7— p-, and JJ = 



1 - ijq/{2 + ^) ' 2q 
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Proof. We proceed analogously to Proposition [T] The dominating singularity of 
(pL^z) is located at z — zi. We use the expansion 

valid for 9 > 0. From this we find 

iPLiziz) ^ CD-\-' /^2'-'rie){i - z)-^ z ^ 1. 

The result now follows from singularity analysis, since the coefhcients of (1 — z)~^ 
asymptotically equal n^^^ /r{d). Once again, an asymptotic expansion to arbitrary 
order can be readily obtained. □ 

Proposition 4. Let N ~ Poisson(A) and A ^ GIG('0,x, §)• Furthermore as- 
sume that X ~ Geo{l,q) with q G (0, 1). Then the computation of the an by the 
recursion (j3.6p is numerically stable. 

Proof. Completely analogous to the proof of Proposition[2l The indicial polynomial 
is w{w — — 2){3w + 2). The root w — leads to a solution whose coefRcients 
grow like ()4.3|) (with 6* = |), whereas the coefRcients of the other solutions grow 
slower. □ 

The approach we have just illustrated works whenever the dominating singulari- 
ties of the differential equation satisfied by fL{z) are regular. In our third example, 
however, the dominating singularity is irregular. In general, it is difficult to say 
anything about the growth order of the power series coefficients of the solutions 
in this case. In our example, though, all solutions of the differential equation for 
(Pl{z) can be expressed in closed form, and their coefficients can be analyzed by 
Cauchy's integral formula and the saddle point method. 

Proposition 5. Let N ^ Poisson(A) and X ^ NBin(i,p) with p £ (0, 1). Then 
the computation of the an by the recursion ()3.8p is numerically stable, and the 
probabilities satisfy 

\l/3„l/6 , 



21/3^3^ 



(1 -p)"n-5/6g^p (^3pl/3(A/2)2/3„l/3 _ 



Proof. We present the proof for p = i and A = 1. The general case yields no 
additional complications. The functions 

±1 



exp ■ 



x/2~z. 

form a fundamental system for the differential equation (|3.7p . Therefore, a funda- 
mental system of the third order recursion p.Sp is given by our a„ , the coefficients 
of exp(— 1/-\/2 — z), and (1,0,0,0, . . .). We have to show that the coefficients a„ of 
(Pl{z) = exp(l/\/2 — z) have the announced asymptotic behavior, and that those of 
exp(— 1/\/2 — z) grow slower. (The additional solution (1, 0, 0, 0, ... ) of (|3.8p can- 
not make the computation unstable, of course.) To do so, we appeal to Cauchy's 
integral formula: 

„ _ 1 /" 'Pl{z)._^ 



2i^./|,|=, 2:"+! 



^ f e-'"'ipL{re')de, < r < 2. 
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We will determine the asymptotics of the integral by the saddle point method [3l [6] . 
To find an approximate saddle point, we equate the derivative of the integrand to 
zero, which leads to the equation 



Clearly, we must have z ^ 2 as n — > oo here. By plugging z — 1 — u with unknown 
u = o(l) into the equation, we obtain u ^ •nr'^l^ . Therefore, we choose the new 
integration contour |z| = r 2 — iit'^I'^ . The dominant part of the integral arises 
near the saddle point, for Q = 0(n^"), where a is a fixed parameter with | < a < |. 
Outside this central part we have 

cos 6* < cos(n-") = 1 - + 0(n"''"), 

and using this estimate in 

\2-z\-^l'^ = (4-4rcos6l + ^2)"l/^ 

we obtain 

< exp(ni/3 _ n5/3-2")(i + o(l)). (4.4) 



exp ■ 



11/2 



To calculate the central part of the integral, we compute the second order approx- 
imation 

(2 - T^'Y^I^ - + ind - \Tt>IH2 + 0(n^/3-3a)^ 

Since 

y" exp(-|n5/3(?2)d0 ~ V?""'^' 
and r~" ^ 2~" exp(in^/'^), we find 

' ^""\--exp(^^)d. ' ""-^^^-^'^^ 



Note that we have shown above that the remaining portion of the integral, where 
ri~" < 16*1 < TT, grows slower, by virtue of the factor exp(— n^/^"^") in (14. 4p . 

Now that we have established the asymptotics of a„ , it remains to show that the 
coefficients, 5„ say, of exp(— 1/-\/2 — z) grow slower. To see this, we use Cauchy's 
integral formula with the same contour as above: 



e-™''exp(-==)d0. 



Here the integrand has no saddle point near z = r, but a bound good enough for 
our purpose can still be deduced. The tail l^j > n"" satisfies the same estimate as 
for (pl{z). Near the real axis, for 6 = o(n""), we have 

I exp(-(2 - re'^)-i/2)| ^ cxpi-n^^^ + f n^/^^^) 

< exp(-ini/3) 

for large n, which shows that the integral over the central part grows slower than 
that for On (it even tends to zero), whence bn = o(a„). □ 
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